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Abstract 

In this paper, we present a systematic variational description to general mass action kinetics of chemical 
reactions with detailed balance by an energetic variational approach. Our approach starts with an energy- 
dissipation law of a chemical reaction system. We show the dynamic of the system is determined by the 
choice of the dissipation. This approach enables us to couple chemical reactions with other effects, such 
as diffusion and drift in an electric field. As an illustration, we apply our approach to a non-equilibrium 
reaction-diffusion system in a specific but canonical set-up. We show the input-output of such a system 
depends on the choice of dissipation via numerical simulation. 


1. Introduction 

Chemistry describes reactions that convert one compound into another. In biological systems, chem¬ 
ical reactions are catalyzed by enzymes and combined to perform many of the functions of life. The 
reactions are combined in networks analogous to networks of electrical circuits: enzymes localize individ¬ 
ual reactions (as well as catalyze them). The products of one reaction move by diffusion (and perhaps 
migration and convection) to become the reactants of another reaction catalyzed by another enzyme in a 
different location. Reactions in biology occur in different physical locations, so the products move as they 
become reactants in a subsequent reaction. 

Here we try to link chemical reactions using the approach of electrical engineering. We describe each 
reaction as a separately defined device (loosely speaking) with an input and output and its own input 
output relations. We do not try to construct a single reaction diffusion system of (for example) partial 
differential equations, just as engineers do not describe an amplifier (for example) by a single grand partial 
differential equation. 

In our approach, the output of one reaction is one of the inputs of another reaction, after the chemical 
species involved diffuse from one location to another. Many types of forces and flows can move the chemical 
species and the partial differential equations of one type of field interact with those of another. We turn to 
the theory of complex fluids, and its variational formulation, energetic variational approaches (EnVarA), 
to describe these interactions because of its success in dealing with otherwise intractable problems. In this 
paper, we aim to embed the classical description of simple chemical reactions into the theory of complex 
fluids, hoping this interdisciplinary interchange will bring new insight, questions, and techniques to both 
fields. 

Since 1950s, there has been a huge amount of work devoted to understand the mathematical structure 
of chemical reaction systems [HU [32i ini eu issi Eli mi Ezi m m osi H 21 usi es us no]- As pointed out 
in J2]: “a biochemical reaction system consists of two parts: (i) a reaction network, and (ii) a choice of 
dynamics.” Depending on the scale of the system, either deterministic or stochastic model are used to 
describe the dynamics of chemical reactions. 
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beisenbe@rush.edu (Bob Eisenberg) 
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The simplest deterministic description to chemical reactions is the law of mass action, which states 
that the rate of a reaction is proportional to the concentrations of the reactants and the proportionality 
constant k is called the rate constant ®H9]. This approach originally arises from the treatment of ideal 
gases oTil . where molecules/atoms only interact when they collide. For example, for a single reversible 
reaction, which occurs in both forward and reverse directions, 

aA + ^B^yC, (1.1) 

k r 


where kf and k r are rate constants for forward and reverse directions, the reaction rate is defined by 


r = rf — r r = fc/[A]“[B ] 7 — fc r [C] 7 , 


( 1 . 2 ) 


according to the law of mass action Q35]. Hence, the rate of changes of [A\, [B\ and [C] satisfy a system of 
ordinary differential equation 

^S. = -a(k f[ Ar [ Bf-k r [CD 

& = ^H{kf[Ar[Bf -k r [CD ( 1 . 3 ) 

^kfiAriBf - kricy) 

At an equilibrium, in which concentrations are not changing, we have 


[^]eij[-®]eg _ t i „ 

[CU ~ V f v ' 


( 1 . 4 ) 


where K eq = is called the equilibrium constant. 

Although the law of mass action has been widely used, as aptly pointed out in Ref. [193, “the law of 
mass action is not a law in the sense that it is inviolable, but rather is a useful model, much like Ohms law 
or Newtons law of cooling. ” Moreover, the law of mass action is mostly used in equilibrium treatments 
of chemical reactions and so does not immediately apply to the biological and electrochemical systems in 
which physical chemistry finds some of its most important applications. Equilibrium treatments assume 
that the system will reach a chemical equilibrium, in which each reversible reaction in the system is in 
detailed balance, that is species have zero flow and zero velocity. 

Nonequilibrium treatments require significant extension to deal with flow, involving new variables and 
the held equations that describe the new variables, along with boundary conditions some within and others 
on the edge of the system. Most of biology and nearly all of electrochemistry occurs in systems with large 
flows and significant dependence on friction. Friction is always present in condensed phases like biological 
solutions, because there is so little empty space between molecules. Innumerable collisions of molecules 
occur in any chemical reaction on the biological time scale. The macroscopic effect of these collisions is 
friction. It has been a challenge to couple the law of mass action with these effects. 

One of the goals of this paper is to provide a seamless and expository extension of the equilibrium 
treatment of chemical reactions to nonequilibrium systems with diffusion and boundary effects. We adapt 
and exploit the theory of complex fluids El GEE] that has dealt with flows in systems with many com¬ 
ponents successfully for many years. These components often each have internal structure and dynamics 
and the theory of complex fluids deals with these as well. 

In our model, the species of ions (and solvent) in electrolyte solutions are the components of the 
complex fluid. We treat all chemically reactive species as ideal gases, so the entropy of the components is 
that of an ideal gas. The nonequilibrium free energy J- of the ideal gas also includes an internal energy 
used to describe their ability to enter into chemical reactions. [For us, a chemical reaction is one in which 
covalent bonds change and the spatial distribution of electrons in the products is not a small perturbation of 
that in the reactants.] Our treatment is always nonequilibrium so J- includes the effects of friction, coupling 
of flows of different species, as well as the effects of multiple fields, e.g., concentration fields (diffusion) 
and electrical fields (migration). The extension to systems with variable temperature is well underway by 
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adopting the methods in [25] ■ We provide a variational version of the theory of complex fluids because 
it guarantees a consistent mathematical formulation, in which all variables satisfy all equations, of all 
fields, and their boundary conditions. Unlike most previous work, which aims to establish the underlying 
energy-dissipation structure for given chemical kinetic equations, our variational approach starts with the 
energy-dissipation law of a chemical reaction system. We show that the chemical reaction are driven by 
the difference of internal energy between “reactants” and “products”, and the dynamic of the system is 
determined by the choice of dissipation, which can be obtained from experimental measurements. 


2. Chemical Reactions in Ideal Gases 

As an illustration, we focus on a single reversible chemical reaction, 

oA + /3B 7 C. (2.1) 

Our framework can be extended to more complicated reaction networks with detailed balance. We assume 
that the entropy of each component in |2.1| can be well described by the ideal gas model without chemical 
reaction [22], so that other types of interactions can be ignored. We are surely aware of the significance 
of those interactions [S] but have to start somewhere. In particular, we follow the usual practice §5) of 
ignoring effects of the shape of the molecule. 

2.1. Free energy and equilibrium, constant 

The key idea is to treat each species in chemical reactions as “ideal gases” with internal energy. The 
internal energy is determined by the atomic structure of the molecule |33] . Then we can adopt energetic 
variational approaches (EnVarA) [23, 3 G33], which are widely used in modeling complex fluids, to show 
that the chemical fluxes, changes of concentration of each component, are functions of the difference of 
internal energy. The internal energy is treated as the internal structure of each chemical species. It is the 
internal energy that might be calculated by quantum chemistry, in the case of no interactions with other 
atoms or molecules. 

The mechanisms of chemical reaction can be very complicated, as the chemical bonds break and form. 
In our approach, we ignore these mechanisms by treating the molecules (which are sometimes bare atoms, 
of course) as objects in a complex fluid that have an internal energy given by their internal structure and 
its internal dynamics. 

These energies can either be measured experimentally, or computed by electronic structure calcula¬ 
tions. The energy of one species of molecules or atoms can be determined by experiments, or in favorable 
cases by calculations of electronic structure, using the techniques of quantum chemical computations. 
Since in a certain sense the internal energy of the electrons of a molecule or atom is electrostatic (once the 
spatial distribution and correlations of the electron orbitals is determined from quantum calculations), it 
is natural to include this term as an energy of a species. Our system then falls naturally in the class of 
complex fluids, studied by the theory of complex fluids [25] , 

What is important here is that we assign an internal energy to each molecule-really to each molecular 
or atomic species-that acts as an input to our further calculations. In our model, chemical reactions 
are driven by the change in the internal energies between reactants and products, along with any forces 
provided by external sources, e.g., boundary conditions. 

In order to simplify the algebra, in the following, we take a = f3 = 7 = 1. Then the reaction becomes 

A + B^C, ( 2 . 2 ) 

We denote the internal energy of A, B and C by C/ 4 , C/g and Uc respectively. So the nonequilibrium free 
energy J- of this system used in the variational theory of complex fluids can be written as 

F{[Al [B], [<7]; Ua, U b , Uc) = j RT([A]{\n[A] - 1) + [B](ln[B] - 1 ) + [C](ln[C] - 1 )) 

+ [A]U a + [B]U b + [C]U c dx, 
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where Ua, Ub and Uc are inputs (parameters). Note that our free energy T is a nonequilibrium quantity 
that depends on friction, coupling between flows, and several kinds of forces. It cannot be derived by 
mathematics alone because it involves a physical model of the effects of friction, of coupling between flows, 
and of the various force fields. We adopt very simple models here to illustrate our approach. Confrontation 
with real experimental data will undoubtedly motivate more complex models. 


The first three terms in (2.3) are free energy of a mixture of ideal gases without chemical reaction, 
which corresponds to the entropy. Indeed, for a mixture of ideal gases with N species, the chemical 
potential of a substance j is expressed by [22, 


Mj = fj° + RT In Xj , 


(2.4) 


where fi° is the reference chemical potential, and Xj is the mole fraction of the substance j. Since the 
chemical potential is defined relative to its value at the arbitrary reference state, we can take /jP = 0 . The 
free energy of the mixture of ideal gases, corresponds to the chemical potential (2.41 with /iq = 0 , is given 

by 


N 


T[x.j] = RT ^Xi(lnXi — l)dx. 


(2.5) 


i=1 


The last three terms in (2.3) are the internal energies stored inside the molecular A, B and C. In the 
case without chemical reaction, since [A], [B] and [C] do not change with respect to time, these terms are 
constants that can be ignored. 

For given free energy ^([A], [B\, [C]; U a ,Ub,Uc) defined in (2.3), the corresponding chemical potential 
of A, B and C are 

VA = -^\=RTln{A} + U A , 

fiB = W] = RT HB] + t/B, (2.6) 

^ c = -^- ] =RT\n[C}+U c . 


At an equilibrium, the chemical potential of both sides of reaction ( 2 . 2 ) are equal, that is 

MA + Mb = Me, 


(2.7) 


so we have 


where A U — Uc — U A 


In 


f [A1 eg [B\ eg \ 
l [C]e, ) 


1 

RT 


{U c -U A - U b ) 


A U 

RT ’ 


( 2 . 8 ) 


Ub is the difference of internal energy between state {A, B} and state {C}. So 


[A] eg [T?] eg 

\C\eq 


A K 


(2.9) 


where K eq is called the equilibrium constant of the reaction (2.2) m- It can be noticed that K eq is an 
exponential representation of the difference in internal (chemical) energy that drives the chemical reaction. 


Remark 2.1. /i^ + — Me known as affinity of a chemical reaction, introduced by De Donder as a 

new state variable of the system. This affinity is the driving force of chemical reaction m in- 


in the literature, the free energy of a chemical reaction system with detailed balance has been written 
down in various equivalent forms pa Husnu si Eg. For the reversible chemical reaction 


k f 

A + B^C, 

k r 


( 2 . 10 ) 
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a conventional form of the free energy (Lapin functional) of the system is given by j.'171 7) 1261 [7] 17] 

T([A],[B],[C}) = - 1) + [R](ln (M) - 1) + [C] In (M) - l)dz, (2.11) 

where (a, 6, c) is the concentrations of A, B and C in a steady-state of system. Another form of free energy 
is given in [17] as 


H\A], [B\, [C]) = J [A]ln(k f [A]) - 1) + [B]ln(k f [B]) - 1) + [C] In (k r [C}) - l)dx, (2.12) 

Although the forms of these free energy seem different, they produce the same equilibrium satisfying 

^=K eq , (2.13) 

eq 

where K eq = fj = ~- The choice of ’best’ form of the free energy depends on experimental measurements 
and other consideration. 



2.2. EnVarA description for chemical reactions 


We are now ready to give an EnVarA description [73] [3] [11] for the kinetics of chemical reactions, 
and to reach towards our goal of showing how the (chemical) internal energy changes the experimentally 
measured variables, namely the number density of products as a function of the concentration of reactants 
and other parameters. 

We start with a brief introduction to an energetic variational approach, which is an extension of a 
variational principle proposed by Rayleigh [158] for purely frictional systems that Onsager tried to extend 
to physical systems in general [221 SO] - During the last decades, EnVarA and its equivalent forms have 
been successfully used to model many complicated systems [73 ini mug. 

The starting point of an EnVarA description is a prescribed energy-dissipation law, which is a conse¬ 
quence of first and second law of thermodynamics. For a closed isothermal system, the energy-dissipation 
law is given by 

^ Jr [z\ =-V[z,z t ], (2.14) 


where z is the state variable, z t is the rate of the state variable, T is the free energy for the nonequilibrium 
system, and T>[z, z t \ is the rate of energy-dissipation. Here we assume the kinetic energy of the system 
can be ignored due to the time scale being considered. The existence of energy-dissipation law (2.14) is a 
consequence of the first law and second law of thermodynamics m- 

For a given energy-dissipation law (2.14), EnVarA provides a general framework to determine the 
dynamics of a complex system through two distinct variational processes: Least Action Principle (LAP), 
taking variation of action functional A = f 0 —Bdx with respect to z, and Maximum Dissipation Principle 
(MDP), taking variation of |V with respect to z t . 

The dissipation T>[z, z t ] is often assumed to take a simple form 0126], 


V[z,z t ] = [Q( y z)z t , z t ), 


(2.15) 


where (.,.) is an inner product, Qz is a positive-definite matrix for given z. The assumption (2.15) 
corresponds to linear response theory in nonequilibrium thermodynamics (75; ;3DI [B] • Then by a standard 
energetic variational procedure, the dynamic of system is given by 


that is 


SB 

(2.16) 

Sz t 6z ’ 

G(r w 

G(z)z t - ^ . 

(2.17) 
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We can extend the standard energetic variational approach to a system, in which the dissipation rate 
T> is not a quadratic function of Zt . Indeed, the dissipation T> can be written as 



V[z,z t ] = (T( 

then 

d / 


&l< 

II 

implying 

V 

II 


Zt 


5JF\z\ 

Sz 

5T[z 


5z ’ 


as a consequence of the chain rule. It can be noticed that if V satisfies (2.15), then 


(2.18) 

(2.19) 


( 2 . 20 ) 


( 2 . 21 ) 


which is consistent with the maximum dissipation principle. 

The main difficulty in applying the above EnVarA to a chemical reaction system is to deal with 
constraints that come from the conservation of chemical elements. For the reversible reaction |2.1[ the 
reaction rate r for this reaction is defined by (see p. 569 in 0 ) 


1 d[A] _ 1 d[H] _ 1 d[C] 

ad t /3 d t 7 dt 

A direct consequence of such definition is that 

^( 7 [A] + a[C])= 0 , A( 7 [£] +/ 3 [C])= 0 , 

if there is no other reaction involving A, B and C. Then 


( 2 . 22 ) 


(2.23) 


7 [A] + a[C\ = Z 0 , 7 [B]+a[C]=Z 1 , 


(2.24) 


where Zq and Z\ are constants determined by the initial amount of A, B and C. Experiments often report 
the concentration (i.e., number density) of A , B and C as a function of time. 

Remark 2.2. For a system that has 3 species, and 1 reaction, we have two conserved quantities. The 
choice of conserved quantity may not be unique in general. For example, for a simple chemical reaction 

H+ + OH - H 2 O, (2.25) 


the conservation of the numbers of H-atoms and O-atoms give us 

[H+] + [OH - ] + 2 [H 2 O] = C H 
[OH - ] + [H 2 O] = Co 


(2.26) 


where Ch cmd Co are numbers of FI atoms and O-atoms. 
is equivalent to the conservation property 


The conservation of (chemical) elements (2.26) 


[H+] + [H 2 O] = Z 0 , [OH - ] + [H 2 O] = Z\ 


(2.27) 


(2.24). 


In order of overcome this difficulty, we introduce a reaction coordinate R(t) to describe a chemical 
reaction system and use the reaction coordinate R(t) as the state variable. The idea of using chemical 
reaction coordinates to describe chemical reactions appeared in 35, 'Tj for both deterministic and stochastic 
descriptions of chemical reactions, but apparently has not yet been used to describe chemical reactions in 
a variational treatment. 
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Roughly speaking, the reaction coordinate accounts for the “number” of each chemical reaction that 
has occurred by time t. We take the system with single chemical reaction (2.2) as an example. By 
introducing the reaction coordinate R(t), the concentrations of A, B and C are given by 


[A(t)] = [A(0)]-i?(t), 
[B(t)] = [B(0)]-R(t), 

[c(t)] = [cm+m, 


(2.28) 


which is a kinematic description of the system that embodies the constraint (2.24). 

By using (2.28), the energy ,F([H], [B], [C];Ua,Ub,Uc ) defined in (2.3) can be written in terms of 
R(t), so the energy-dissipation law can be reformulated as 


-F[R- U A , U b , U c } = -T>[R,Rt], 


(2.29) 


with a proper (but not unique) choice for the dissipation V[R,R t ]. Here, the energy T[R;Ua,Ub,Uc] 
is a functional of R(t) with Ua, Ub and Uc as inputs (parameters). The energy dissipation law (2.29) 
embodies the constraints 12.241 

We can derive the equation of R(t) by the generalized energetic variational approach 2.20 Note by 
the kinematic equation (2.28), we have 


d[A] 


— —Rt, 


d [B] 


= -Rt 


d[C] 


= Rt 


d t ’ dt d t 

so the reaction rate law r is uniquely determined by the choice of the dissipation T>[R, Rt]- 


(2.30) 


2.2.1. Law of mass action 

In the law of mass action, the reaction rate is given by 

r = k f [A][B\ - k r [C] 


(2.31) 


where kf and k r are rate constants for the forward and reverse directions, the corresponding equilibrium 
constant is 


Keq = 


eq [-^j eq k r 


\C] 


eq 


V 


(2.32) 


as 


The law of mass action can be derived from the energy-dissipation law (2.29) by taking the dissipation 

Rt 


V[R,R t } = RT Rt In 


k r [C\ 


+ 1 


(2.33) 


where k r is the reverse rate constant of this reaction. Indeed, the energetic variational procedure gives 


RTln 


Rt 


k r [B] 


+ 1) = -j^T[R-,Ua,U b ,U c }. 


Note 


SR 


R[R; Ua, U b ,Uc ] = (—- Ms + Me) 


= RTln 


(xm\)- UA ~ UB+Uc 


by using (2.6). So we have 


In 


Rt 


k r [B] 


+ 1 


=/m - ^ 


V [c] ) 


RT 


(2.34) 


(2.35) 


(2.36) 


where the right-hand side is determined by the difference of internal energy A U between state { A , B} and 
state {C}. 
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Although (2.36) looks complicated, direct computation gives 


r = R t = k r [C] - l) = k f [A][B] k r [C ], 


where 


is used. 


A U k r 

A - = eRT =v f 


(2.37) 

(2.38) 


Other forms of the dissipation functional are possible. As a generalization of (2.33), we can consider 
a more general form of the dissipation 


V[R,Rt] = r n (R)d t R\n(r l2 {R)R t + 1), 


(2.39) 


where 771 (i?) > 0 and 772 (i?) > 0, then V(d t R) > 0 for the admissible R t . In this case, by (2.20), we have 


vi(R) Hv 2 (R)R t + 1 ) = u a, u b ,u c ] 


= RTln 




(2.40) 


In this case, the forward and reverse rate kf and k r can depend on the concentrations of [A], [B\ and [C]. 
For instance, if 771 ( 1 ?) = RT and 772 ( 1 ?) = 1, then 


where 


r = Rt = ' 1 = kf([C])[A)[B] k r ([C])[C], 


fc/ ([ C] )- [C\K eq ’ fcr([<7]) ~ [C\ 


(2.41) 

(2.42) 


2.2.2. Linear Response Theory 

In nonequilibrium thermodynamics [291 1301 [6] , it is often assumed that the dissipation of total energy 
is a quadratic function of “rate” of state variables, which is known as linear response theory. 

It is important to note the distinction between the linear response used so widely in chemical kinetics 
and the linear response used in engineering. In engineering, an operating point in a complex multidimen¬ 
sional space describes the nonlinear properties of a system. This point is usually remote from, and has little 
to do with the properties of the system with no flows, i.e., no sources of energy, mass, charge or current, 
that is called equilibrium in chemistry. Linearization around that operating point is an essential tool in 
understanding most engineering systems, even those as nonlinear as digital modules (say shift register to 
be specific) in a digital computer. But systems without flows are rarely considered in engineering since 
they have so little to do with the systems that actually perform engineering function. Think of an amplifier 
without a power supply, or an automobile engine with water in the gas tank. 

In our case, the linear response theory of chemistry gives us a form of the dissipation term 

V[R,R t ] = V (R)\R t \ 2 . (2.43) 

Then a standard variational procedure gives us 

v{R)Rt = ~muA,u B ,Uc] 

(Jcl\ 

U MR}) 
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= RTln 


-A U. 


(2.44) 










By choosing r](R) = RT, the reaction rate is given by 


r = Rt = In 


= In 


[C] \ A U 
[A][B]J RT 
1 [A][B]\ 
K eq [C] )' 


(2.45) 


which is much more complicated than the classical reaction rate (2.311. The different reaction rates 
obviously will predict different dependence of rate on concentration and different time courses of the 
chemical reaction. We reiterate that the choice of dissipation function is a scientific not mathematical 
decision. The dissipation function must be chosen to fit data and/or to describe a model or simulation, 
as well as the mathematical requirements of EnVarA. The linear response theory arise from the near 
equilibrium assumption. For the chemical reaction, we have Rt ~ 0 near the equilibrium, then the Taylor 
expansion gives us 

m{R), 


Vi(R) In (ri2{R)Rt +1) 


m(R) 


I Rt 


(2.46) 


So one can view the dissipation (2.43) as a linear approximation near equilibrium to (2.39). 

Remark 2.3. The law of mass action gives a simple form of the reaction rate r in terms of [A] and [B], 
however, the dissipation in terms of R and Rt would become complicated (See equation 2.33). On the other 
hand, if the dissipation is taken to be simple that described by linear response theory, the the reaction rate 
r becomes complicated (See eq. (3.12) ). 


Remark 2.4. Note K eq is defined as an exponential representation of the difference in internal (‘chemical’) 
energy, measurements of just K eq cannot distinguish these dissipation mechanism. Measurements of the 
concentration dependencies of the forward and backwards rates can distinguish them. However, it is difficult 
to measure the forward and backwards rate respectively. 


3. Diffusion and boundary effects 

In our EnVarA description, the dynamic of a chemical reaction is determined by the choice of dis¬ 
sipation. The dissipation is best specified in the context a specific experimental set up. Such setups are 
chosen to give reproducible input-output functions that are useful in applications. 

Our analysis will show how the setup determines a dissipation function in a particular case, as an 
example of how that might be done more generally. The setup of the system is shown in Fig. |3.1[ 





channel 





Figure 3.1: Setup of system 


In this system, a narrow channel connects two bath. In left bath, the average concentration of A, B 
and C are maintained by the boundary condition. We only consider the chemical reaction 

A + B^=^C (3.1) 

happens inside the channel. The species in the left bath are sources. We are concerned with the change 
of amount of C in the right bath, and assume both A and B cannot diffuse into the right bath. The flux 
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of C or the rate of change of amount of C in the right bath is our output. The concentration in the left 
bath are our inputs. The chemical reaction is the “transfer function”. 

Mathematically, this system can be modeled by coupling the reaction with diffusion in the framework 
of energetic variational approach. Since the channel is very narrow, we can treat this problem as an 
onedimension. Let c A ,c B and cc be concentrations in the channel [—e, e]. Due to the above setup, we 
have boundary conditions for 

ca{~c, t) = c\, d x c A (e , t) = 0 , 

c B (~e, t) = 4, d x c B {e , t) = 0, (3.2) 

d x c c {-e , t) = 0, c c (e, t) = c b c . 

Here we impose the non-flux boundary condition for cc on the left-end of the channel, and Dirichlet 
boundary condition on the right-end of the channel. The boundary conditions of c A and c B are Dirichlet 

boundary condition, which are inputs of our system. Since cc satisfies the Dirichlet boundary condition on 

the right-end of the channel;, the amount of C diffuse into right bath by time T can be roughly computed 
by 

R(x, T)dx — J (c(x,T) — co(x))dx, (3-3) 

which is the output of our system. The flux of C or the rate of change of amount of C in the right bath 

can be computed as —C out . 

at 

We choose reaction coordinate R and fluxes J Al J B and Jc as state variables of this system. The 
kinematic equations for c Al c B and cc are given by 

c A (x,t) + V • J A {x,t) = c° A (x) - R(x,t), 

c B (x, t) + V • J B (x, t) = Cg(x) - R(x, t), (3.4) 

cc(x, t) + V • Jc(x, t) = Cq(x) + R(x , t), 


C out (T) = J 


where the initial conditions c A (x), cP rj (x). and c° G {x) satisfy the boundary conditions. The energy- 
dissipation law of the system can be formulated as 


d 

df 


J c A (In c A - 1 ) + c B (In c B - 1 ) + c c (lnc c - 1 ) + c A U A + c B U B + c c U c dx 
- j D(R(x,t),d t R(x,t)) + VA\d t J A \ 2 + r] B \d t J B \ 2 + r) B \d t Jc\ 2 dx, 


(3.5) 


where the energy part is exactly the same as the energy in the previous section with U A , U B and Uc as 
inputs, the dissipation part includes dissipation with respect to dtR and dtJ■ 

By using the generalized energetic variational approach |2.20| we have 


1 


D(R, d t R ) = fJ. A + f’B — fJ-c 


(3.6) 


r)i(J)d t Ji = -Vni(R,J), i = A, B, C, 


where fit (i = A, B , C) is the chemcial potential for each species, as given in (2.6). By choosing 7) l = A(i = 
A,B,C) and combining with (3.4), we can obtain a standard reaction-diffusion equation 


( d t c A = V • (Vc A ) - r(x, t ) 

< d t c B = V • (Vc B ) - r(x,t) (3.7) 

[ d t c a = V • (Vc c ) - r(x,t), 


where r(x,t) is determined by the choice of D(R,dtR ) as discussed in subsection 2 . 2.1 and 2 . 2.2 
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3.1. Simulations and discussion 

In this subsection, we show some simulation results for the above system with different choices of 
dissipation. We fix e = 0.1 throughout this subsection. As mentioned previously, the inputs of the system 
are concentrations of A, B and C in the left bath, which gives us the Dirichlet boundary conditions of A , 
B and C in the left-end of the channel. We assume that 


C A ~ C B — C 0 


(3.8) 


such that Co is the single input of our system. The initial concentrations of A , B and C in the channel 
are constants c° A (x) = c° B {x) = Co and c° c (x) = 0.1. The output of system is C out defined in 3.3 which 
roughly correspond to the amount of C diffuse into the right bath, or , the approximated flux to the 
right bath. 

We fix K eq = 0.1 through this section and show the input-output relation depends on the choice of 
dissipation. We focus on two types of dissipations, as discussed in subsection 2.2.1 and 2.2.2 The first 
form of dissipation is taken as 


D\ (R, d t R) = rn(R)R t ln(r] 2 (R)R t + 1), 


(3.9) 


which corresponds to generalized law of mass action. The second form of dissipation is taken as 


D 2 (R : d t R)=rj 0 (R)\Rt\ 2 , 


(3.10) 


which corresponds to linear response theory. We take rjo(R) = r]i(R) = ?y 2 (R) = 1 in our simulations. In 
general, r)i(R) is also a changeable variable that depends on R. 


Fig. 3.2 shows the output C out (t) as a function of input cq when t = 1 for the two choices of dissipation. 
It can be noticed that for small c 0l the output are close between two choices for dissipations. However, 
the output for dissipation (3.91 is much larger than that for the dissipation (3.101 when cq is large. 



Figure 3.2: The output CoutR) as a function of input co when t = 1 for dissipation { | 3 .9 [ ) [circle] and dissipation ( |.3. 1 0[ i 
[square]. 


Formally, from the computations in subsection |2.2.1| and |2.2.2[ we know 


r = Rt = In 


( 1 IMB}\ 
\K eq [C] ) 


(3.11) 
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for the dissipation (3.9), while 


r = R t = 


1 [MB] 

K eq [C] 


- 1 


(3.12) 


for the dissipation (3.10). For Co = 0.1, the system is at the equilibrium, so C out = 0. When Co is large, 
the dissipation (3.101 determines large reaction rate. 

We also consider the 4 C ou t as a function of t for two choices of dissipation for various of Co- The 
results are shown in Fig. |3.3[ The time courses in Fig. |3.3| show that for both dissipations and different 




Figure 3.3: A Cout as a function of t for two choices of dissipation for various of co (co = 1, 0.5 and 0.25 from top to bottom 
in each figure), (a) Dissipation (|3.9[), (b) Dissipation (|3. 10 [i. 


inputs, jf Cout tends to a constant, which is a function of input for a given dissipation. 

The above simulation results suggest that one can determined dissipation through experimental mea¬ 
surements and solving the inverse problem. 


4. Summary 

In this paper, we establish a variational structure to a chemical reaction system with detailed balance, 
which enable us to couple chemical reactions with diffusion, boundary effects and other effects of systems. 
In our approach, the dynamics of a chemical reaction is determined by the choice of the dissipation. The 
classical law of mass action can be derived via a particular form of dissipation. 

As an illustration, we apply our approach to a nonequilibrium system with boundary effects and show 
that the input-output of such a system depends on the choice of the dissipation via numerical simulation. 
Since such a system is experimentally reliable, it indicates that the dissipation functional can be obtained 
by experimental measurements and studying an inverse problem. 

There are potential applications in our energetic variational description to a system involving chem¬ 
ical reaction. For instance, new numerical schemes can be designed based on our energetic variational 
formulation. It is often a challenge to preserve the non-negative and conservation properties mm- Since 
we have embodied the conservation properties by choosing the proper state variable, such conservation 
properties can be naturally preserved for the numerical methods based on our energetic variational form. 
Moreover, such energetic variational formulas also enable us to design Lagrangian-Eulerian schemes for 
reaction-diffusion systems by applying some recently developed methods for general diffusions [I8ll4l 24ll3]. 
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